//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
/*************************************************************************
 * template<class T>
 * class BoostLapackSymSolver
 *
 * Input Matrices = Lower Triangle Packed (in fortran)
 * using triangle_matrix<T,upper> 
 *
 * Expert driver to calculate M eigenstates, M < N 
 *************************************************************************/
#ifndef OHMMS_OHMMSMATRIX_ESSL_SYMSOLVER_PACKED_H
#define OHMMS_OHMMSMATRIX_ESSL_SYMSOLVER_PACKED_H

#ifndef EIGENSPACEDIMENSION_H
#include "EigSolver/EigenSpaceDimension.h"
#endif

extern "C" {
  void dspev(int*, double*, double*, double*, int*, int*, double*, int*);
  void zhpev(int*, complex<double>*, double*, complex<double>*, int*, int*, double*, int*);

  void sspev(int*, float*, float*, float*, int*, int*, float*, int*);
  void chpev(int*, complex<float>*, float*, complex<float>*, int*, int*, float*, int*);

  void dspsv(int*, double* , double*, double*, int* ,int* ,int *, double* , int*);
  void dsygv(int*, double*, int*, double*, int*, double*, double*,
             int*, int*, double*, int*);

  void dgegv(int*, double*, int*, double*, int*, double*, double*, double*,
            int*, int*, double*, int*);
}

template<class T, class Ham_t, class EigVec_t>
class EsslSymSolver {

public:

  typedef typename OHMMS::type_traits<T>::value_type value_type;
  typedef typename OHMMS::type_traits<T>::real_type  real_type;
  typedef Ham_t    HamiltonianMatrix_t;
  typedef EigVec_t EigenStates_t;
  typedef vector<value_type> WorkArray_t;
  typedef vector<real_type>  RealArray_t;
  typedef vector<int>        IndexArray_t;
  typedef typename WorkArray_t::size_type   size_type;
  EigenSpaceDimension& States;       
  
  EsslSymSolver(EigenSpaceDimension& dim, bool nonortho=false): 
  States(dim){ 
    LOGMSG("Creating ESSL SymSolver")
  }


  ~EsslSymSolver() { }

  inline real_type eigval(size_type i) const { return eigVal[i];}

  char Jobz;
  int iopt,naux,N;
  RealArray_t workd; // array of 3N+1
  RealArray_t eigVal;// eigen values
  RealArray_t eigValI;// eigen values

  inline bool solve(double*  a, double *z) {
    int M = States.M;    
    check_workspace(States.N, 3, 0);
    if(N/M > 5) {
      dspsv(&iopt, a, &eigVal[0], z, &N, &N, &M, &workd[0], &naux);
    } else {
      dspev(&iopt, a, &eigVal[0], z, &N, &N, &workd[0], &naux);
    }
    return true;
  }

  inline bool solve(double* ap, double *bp, double* z) {
    check_workspace(States.N, 3, 0);
    dsygv(&iopt, ap, &N, bp, &N, &eigVal[0], z, &N, &N, &workd[0], &naux);
    return true;
  }

  inline bool solve(complex<double>* a, complex<double> *z) {

    check_workspace(States.N, 6, 0);
    zhpev(&iopt, a, &eigVal[0], z, &N, &N, &workd[0], &naux);
    return true;
  }

  inline bool solve(complex<double>* ap, complex<double>* bp, complex<double>* z) {

    return false;
  }

  void check_workspace(int n, int workd_size, int rworkd_size) {
    N = n;
    naux = n*workd_size;
    if(workd.size() < naux) workd.resize(naux);
    if(eigVal.size() < n)   eigVal.resize(n);
  }
};

#endif

/***************************************************************************
 * $RCSfile: EsslSymSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: EsslSymSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

